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In this work we consider vortex lattices in rotating Bose-Einstein Condensates composed of two 
species of bosons having different masses. Previously [l] it was claimed that the vortices of the two 
species form bound pairs and the two vortex lattices lock. Remarkably, the two condensates and 
the external drive all rotate at different speeds due to the disparity of the masses of the constituent 
bosons. In this paper we study the system by solving the full two-component Gross-Pitaevskii 
equations numerically. Using this approach we verify the stability of the putative locked state which 
is found to exist within a disk centered on the axis of rotation and which depends on the mass 
ratio of the two bosons. We also show that an analytic estimate of this locking radius based on a 
two-body force calculation agrees well with the numerical results. 



I. INTRODUCTION 

One of the most striking manifestations of the 
quantum- mechanical nature of superfluids under rotation 
is the formation of vortices Q, Q • Perhaps the most nat- 
ural arena to controllably study the physics of vortices 
are Bose-Einstein Condensates (BECs) of alkali atoms 
0, S Q • For the simplest case where the condensate is 
composed of a single type of atom without spin degrees 
of freedom, a triangular lattice is formed 0, Hf • On the 
other hand, for multicomponent systems (composed of 
mixtures of atoms or spinor condensates, for instance), 
the order parameter has additional degrees of freedom 
resulting in more complex vortex lattice structures (see, 
for instance, 0, EO, El E3, E3 ) • 

For the classic problem of an ideal fluid in a con- 
tainer rotating at rate fi, the steady-state local velocity 
is v = ft x r where r is the distance from axis of rota- 
tion. Since this velocity has everywhere a nonvanishing 
curl (|V x v| = it is not a permissible flow for a su- 
perfluid which is supposed to be inherently irrotational 
(v = W with the phase of the SF order parameter). 
However, superfluids are well-known to mimic the clas- 
sical rigid-body rotation on average by forming a vortex 
lattice where the density of these vortices is given by the 
Feynman relation 0, Q 
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where m is the mass of the constituent bosons and Q is 
the rate the superfluid is rotating which is equal to the 
rotational rate of the walls of the container. 

In a previous work [1], it was considered how the situ- 
ation described above generalizes to the problem of two- 
component BECs composed of atoms having different 
masses. More specifically, Eq. (pQ) naturally generalizes 
for two-component systems to 

raiQi 2 ^2^2 
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where mi and vrt2 are the masses of the bosons in the two 
constituent condensates and Vt\ and ^2 are angular rates 
at which the two superfluids are rotating. For the case 
where there is a negative interspecies scattering length, 
the attraction between species will lead to an attractive 
interaction between vortices of the two species. When 
this interaction is sufficiently large one has the situation 
where the vortices form bound pairs, forcing the densi- 
ties of the two vortex lattices to be the same: p\ « p\. 
For this case, Eqns. ([2]) imply that the two superfluids 
(taking without loss of generality m\ > 777,2) will rotate 
at different speeds Qi < 0,2- This counterintuitive state 
results from the quantum mechanical nature of the su- 
perfluid and has no analog in the classical fluid case. 

In [l| the existence of this state was argued by mak- 
ing an ansatz for the short-ranged interspecies interac- 
tion and performing a two-body force calculation using it. 
This gave a quantitative prediction for the distance from 
the center of the condensate at which the vortex pairs 
become unbound, resulting from the growth of the Mag- 
nus force, which is referred to as the locking radius. The 
goal of the current paper is to test these arguments by 
numerical integration of the full two-component Gross- 
Pitaevskii equation. We will verify that such locked 
states are stable for a range of parameters. Furthermore, 
we will see that the analytic prediction for the locking 
radius agrees well with the numerical results. 

This paper is organized as follows. First in Sec. Ill Al 
we provide definitions and set the notation for the treat- 
ment of the Gross Pitaevskii equation. In Sec. Ill Bl we 
summarize the derivation of the locking radius previously 
given in [l| which is based on a two-body force calcula- 
tion. Then in Sec. [Till we describe the split-operator tech- 
nique utilized to propagate the Gross Pitaevskii equa- 
tions in imaginary time. The main results of the paper 
are presented in Sec. [iVl Here we provide the vortex 
lattice structures determined numerically, and compare 
them with the estimate for the locking radius. Finally, in 
Sec. |V]we provide a discussion of potential experiments 
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to realize this effect and then conclude. 
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II. BACKGROUND 

A. Two-component Gross-Pitaevskii Equations 

Our analysis starts with the two-component Gross 
Pitaevskii energy functional in the frame of reference 
rotating at angular rate Q^. This is given by E = 
Ex+ E 2 + E 12 where 

£1 = j d 2 r A-|V^i| 2 + Vm! + X -g x n\ - Q d ^L z ^ , 

(3) 
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and 
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In these equations, V\ and V2 are the confining potential 
of the BECs which we will take to be harmonic. The 
intraspecies and interspecies scattering strengths are de- 
fined as #i 5 2 and g\2 respectively. The angular momen- 
tum operator, as usual, is defined as L z = xp y — yp x 
where p XiV = —iftd x ^ y . 

Varying this energy with respect to ^1 and ?p2 and 
introducing the chemical potentials fi\ and fi 2 to enforce 
particle number conservation gives the two-component 
Gross-Pitaevskii equations 

Ti 2 o 
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M2^2 = V 2 ^2 + ^2^2 + g2U2^2 + #12^1^2 ~ ^£^2- 
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In Sec. we will describe how these coupled equations 
are solved numerically to find minima of the energy E. 



B. Inter-species vortex attraction and locking 

In this section, for completeness, we provide an es- 
timate of the locking radius of the vortex-bound state 
based on a two-vortex calculation. This calculation is 
presented in more detail in [![. We first consider the 
state for the case where the interspecies interaction is 
large and all of the vortices of one species are bound with 
that of the other due to the strong short-ranged attrac- 
tive force. A bound pair of vortices is depicted in Fig. [H 
The vortex binding causes the two superfluids to rotate 
at different rates (because of the different masses and the 



FIG. 1: A bound pair of vortices occurring in the locked state 
composed of a mixture of BECs with m\ > mi. The vortex 
in the heavier species is denoted with an 'x' while that in 
the lighter species is denoted with an 'o'. The center of the 
condensate is taken to be to the left of this bound pair. The 
attractive short-ranged interspecies force F^ r serves to bind 
the vortex pairs together. This is counterbalanced by the 
Magnus force F^ag which increases from the center of the 
condensate. Note that the healing length of the superfluid is 
larger than the sphere representing vortices in this figure 



Feynman relation), creating a Magnus force which tries 
to rip the bound pair apart. The Magnus force is bal- 
anced by the short-ranged interspecies vortex interaction 
resulting from the overlap of the vortex cores. However, 
since the Magnus force grows linearly with the distance 
from the center of the condensate, it will eventually over- 
come the short-ranged interspecies attraction. The point 
at which this occurs we refer to as the locking radius. 

We will now put the previous arguments on more quan- 
titative footing. A vortex sitting at rest in a superfluid 
flowing at velocity v will experience a force perpendicular 
to the flow 



Fmag = 27rftn v x ft, 



(8) 



the so-called Magnus force [8|], where k is a unit vector 
centered on the vortex pointing out the plane. Assuming 
the system is composed of entirely locked vortex lattices, 
we have for the vortex densities = pv J > . This, via the 
Feynman relations, gives 



(9) 



where Qi and ^2 are the angular rotational rates of the 
two superfluids. Since mi > m2 the superfluids will ro- 
tate at different rates which will lead to the Magnus forces 
pulling the bound pair apart. 

The short-ranged interspecies interaction arising from 
£?i2 defined in Eq. ([5|) will counteract the Magnus force. 
We refer to this as the "restoring force". In order to 
obtain an analytic expression for this interaction, we take 
a Gaussian ansatz for the density profile about a vortex. 
Specifically, for a vortex in species a centered at ro, we 
take 



n a (r) = n£(l 
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where \ a is on the order of the superfluid coherence 
length [l|]. We take 



A a = 1.781& 



(11) 
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for the value of this length parameter. This is a slight 
modification of the analysis in [l[ where the simpler case 
of X a = £ a was taken. We choose the the nonuniversal 
constant in Eq. (jTTJ) so that the ansatz in Eq. ([T0|) pro- 
vides a better fit to the density surrounding a vortex. For 
a more detailed discussion of this, see the Appendix D 

Inserting this ansatz for two vortices separated by dis- 
tance d into Eq. ([5]) we find 
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where we have dropped terms which do not depend on the 
vortex separation. The interspecies force immediately 
follows from the derivative of this interaction energy and 
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Balancing the forces on each vortex in the frame of refer 
ence rotating at the drive frequency, we have i^ag 
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and i^ ag = F r 2 str , as is illustrated in Fig. [TJ Since the 



restoring force acting on either species has the same mag- 
nitude we have that F^ = F^ . The Magnus force on 
species a in the frame rotating with the vortex lattice at 
frequency Q v is given by 



F« ag = 2nhn«\{l a -n v \r 



(14) 



which grows linearly with the distance from the center of 
the condensate r. Also, note that the restoring force will 
not depend on the position in the condensate. A bound 
vortex pair will become unstable when the Magnus force 
is equal to the maximum possible value of the restoring 
force. This can be worked out to be 
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A more detailed derivation of this expression can be 
found in Ref. Q. 



III. NUMERICAL METHODS 

Since there are only rare occasions when the time- 
dependent Gross-Pitaevskii equation permits analytic so- 
lutions, numerical simulation is often the method of 
choice for theoretically studying Bose-Einstein conden- 
sates (for a recent account numerical solution of the GPE, 
see Ref. [14|). In this section, for simplicity, we will only 
consider the single-component case, noting that the gen- 
eralization to the two-component case is straightforward. 
To this end, the equation we wish to solve is 
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which describes the evolution of ip in imaginary time, 
r = it. Under long enough evolution ip will relax to the 
ground state of the Gross-Pitaevskii energy functional. 
In the above equation H is given by 
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We use a split-operator method to evolve the order 
parameter ip as in Eq. (fT6|) . The idea behind the split- 
operator method is to approximate the evolution op- 
erator through imaginary time interval Ar, U(Ar) = 
e~ HAr , by a product of terms which are easily diagonal- 
izable. Neglecting for the moment the rotational term 
in Eq. (fTTj) . H can be written as the sum of two terms, 
H = T + V, where T = -£v 2 and V = V tlap + <#| 2 . 
These terms are easily diagonalized in momentum and 
position space respectively. The wave function ip can 
then be advanced in time by Ar by 



i/j(t + Ar) 
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which is accurate to second order in Ar. The order pa- 
rameter can then be evolved by taking successive Fourier 
(and inverse Fourier) transforms of ip and multiplying 



by the factors e~ ^ r , e , and e ^ respectively. 
Such Fourier transforms account for the bulk of the com- 
putational cost in this algorithm, thus using the efficient 
Fast Fourier Transform algorithm is crucial. 

A complication in the above occurs due to the nonlin- 
earity of the GPE. That is, V in our above prescription 
for time evolution depends on the density n = |^| 2 , and 
it is at first unclear for what time this quantity should 
be evaluated. It is shown in [TEj] that, provided we use 
the most updated version of the time-dependent density 
n = |^| 2 , Eq. (fT8|) will retain its second order accuracy. 
The final complication occurs from the rotational term 
in H which is 
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We have neglected this term thus far since it is diag- 
onalized in neither position nor momentum space and 
therefore cannot be included in either T or V. However, 
we note that R commutes with both T and V so we can 
write 



i/j(t + Ar) 
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e -±TAT e -RAT e -VAr e -±TAT 



V(r). 



Then we can perform a similar split-operator decompo- 
sition of the additional term as 



^hVLxpy At p — hflypx Ar ^hflxp y At 



(21) 



Evolution of ip by this factor can then be performed by 
taking the partial Fourier transform of ip, that is trans- 
forming over the x variables but leaving the y variables 



unchanged (or vice- versa). This completes the overview 
of the numerical method used to solve the GPE. As stated 
before, the generalization to the two-component case is 
straightforward. 



IV. RESULTS 



density) one finds 



Rtf 



(23) 



where r[? is Eq. (fT5|) evaluated for parameters at the cen- 
ter of the trap. This equation can then be solved for r c 
to obtain the renormalized value of the locking radius 



Next, we discuss the results of the numerical simula- 
tion. We first provide the parameters which were used 
for the computations. For simplicity, we restrict our at- 
tention to the simplest case where g\ = g2 , and we fix the 
interspecies scattering strength such that \g\2\/gi = 2/3. 
Furthermore, we take the number of particles in each 
species to be the same: N\ = N2. We set the dimension- 
less parameter defined as g = ^g±Ni to be g = 2 x 10 4 
which is in line with values from typical experiments [2l| . 
We take the two trapping potentials to be harmonic and 
adjust their curvatures ^1,^2 so that the density profiles 
of the two species have the same Thomas- Fermi profiles. 
Finally we rotate the system at 0.9 times the critical rate 
at which the condensate becomes unstable due to cen- 
trifugal forces. We discretize the system on a 200 x 200 
grid, and propagate the system in imaginary time inter- 
vals of At = 0.01^. 

We first consider the simplest case where the masses 
of the two species are the same. For this state we take 
the initial wavefunction to be a perfect triangular lat- 
tice of vortices with density given by the Feynman rela- 
tion, Eq. (pQ). This structure is then relaxed by evolving 
the wavefuntions in imaginary time using the methods 
described in Sec. HTT1 As expected these vortex lattices 
remain fully locked. The relaxed structures show small 
deviations from the perfect triangular initial structure 
due to the effects of the trap [l6| . The density profiles 
of these are shown Fig. [2] in panels (la) and (lb). The 
positions of the vortices are determined by analyzing the 
phases of the relaxed wave functions. Using this relaxed 
structure as the initial state, we change the mass ratios 
and propagate the wavefunctions in imaginary time un- 
til convergence. Specifically, we consider the ratios of 
m\/m2 = 1.2, 1.4. and 1.6 as shown in Fig. O 

To compare these numerical results to our estimate 
for the locking radius described in Sec. IIIB[ we need to 
tailor Eq. (|T5j) to the case of a harmonic trap. We take 
the density profiles used in Eq. (jT5j) to have the Thomas- 
Fermi form: 



nio = n 1 
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(22) 



where no is the density at the center of the trap and 
Rtf is the Thomas-Fermi radius (note that we are only 
considering the case when the two condensates have the 
same radius). Inserting this profile into Eq. ([T5]) (taking 
the correct dependence of the coherence lengths on the 
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This indicates that near the center of the condensate we 
will have r c « as expected. Also, the locking ra- 
dius will never exceed the radius of the condensate as 
expected. The reduction of the bare value of the locking 
radius can be qualitatively understood as follows. The 
Magnus force is proportional to the superfluid density 
while the restoring force is proportional to this density 
squared. Therefore near the edge of the condensate where 
the density is considerably smaller than its value at the 
center, the Magnus force will be favored thereby sup- 
pressing the locking radius. 

Shown in the third column of Fig. [2] are the positions of 
the vortices of the two species, labeled with x's and o's. 
The width of these labels are roughly the size of the co- 
herence length of the condensates. Superimposed on this 
is the locking radius predicted by Eq. (|24|) (using Eq. ([24]) 
for the bar locking radius) shown as a dotted line. This 
shows that the analytic results provide an excellent esti- 
mate of the locking radius. Note that due to the strong 
interactions between the two condensates, an unbound 
vortex in one species will create a local minimum in the 
other. Such features can be seen in columns c and d of 
Fig. [2j These local depletions should not be mistaken for 
vortices which are defined by the phase behavior of the 
wavefunctions. 



V. EXPERIMENTAL CONSIDERATIONS AND 
CONCLUDING REMARKS 

The main requirement for realizing the locked state is 
a BEC composed of a binary mixture of atoms having 
different masses and a negative scattering length. Such a 
transition could be tuned with an interspecies Feshbach 
resonance which have been found in Li-Na [TtJ and Rb-K 
[l8[ mixtures. Such mixtures have respective mass ratios 
of 3.3 and 2.2. Another promising experimental system 
are mixtures of two isotopes of a particular atom. For 
instance, the interspecies scattering lengths of different 
species of Yb have been analyzed in [19| and are often 
found to be negative. Since the mass ratios for different 
isotopes are closer to unity, having a strong attractive 
interaction (often requiring a Feshbach resonance) is un- 
necessary to reach the vortex locked state for this case. 

We also note that these results are closely related to the 
experiment described in (2(j. Here a single-component 



FIG. 2: Relaxed vortex lattices for several mass ratios. Rows 1 and 2 are the density profiles of species 1 and 2 respectively 
as a function of position. Superimposed over these images are the vortex positions marked with red x's and black o's. Row 3 
gives the positions of both vortex species for comparison. In this row, the dotted circle is the estimate for the locking radius 
based on the two-body calculation showing good agreement with the numerics. Columns a, b, c, and d are for mass ratios of 
mi/rri2 — 1.0, 1.2, 1.4, and 1.6 respectively. 



BEC is stirred by a rotating optical lattice which acts as 
vortex pinning sites. When the optical lattice is rotated 
at the speed for which the density of the pinning sites 
matches the density of the vortex lattice predicted by 
Eq. (pQ), a completely locked state is observed. Away 
from this resonance, a similar analysis to the above will 
predict a disk of bound vortices. 

In conclusion, we report the confirmation of the puta- 
tive vortex locked state proposed in [l|. For this state, 
the two superfluids and the stirring potential all rotate 
at different rates, exhibiting an unusual effect due to the 
quantum mechanical nature of superfluids. In this pa- 
per, we showed that such a state exists within a disk 
centered on the axis of rotation and whose size agrees 
well with an analytic estimate. Note that our numeri- 
cal analysis did not assume anything about the vortex- 
vortex attraction (unlike our theoretical analysis, which 
assumes Eq. ([T2]h and evolves the Gross-Pitaevskii equa- 



tions directly). Our results (both analytical and numer- 
ical) rely on approaching this state from the fully locked 
state. Experimentally, this is probably most easily real- 
ized by controllably adjusting an interspecies Feshbach 
resonance. Alternatively, one can use an optical lattice 
to control the effective masses of the atoms by varying 
the lattice depth. 



Acknowledgments 

We would like to thank M. Porter and H.-P. Biichler 
for collaborations on related previous work. We would 
also like to thank L. Baksmaty for valuable advice on 
numerical methods. This work was supported by the 
Sherman Fairchild Foundation (RB); the Caltech SURF 
program (EC); and the Packard and Sloan Foundations, 
the Institute for Quantum Information under NSF grants 



6 




The density n = f 2 resulting from the numerical solution 
of this equation is shown in Fig. [3j 

The numerical solution shows that the density be- 
havior close to the vortex core (r <C £) is n(r) « 

0.340 no (j^J . On the other hand, the far distance be- 
havior is found to be n(r) — no ^1 — (j^j ^. To make 

our work amenable to analytic treatment, we take the 
following ansatz for the vortex profile: 

n(r) =n (l-e- r2 / A2 ) (A.2) 



FIG. 3: Solid line: the density profile for a single vortex at 
r = found from numerically solving Eq. (|A.1[) . Dashed line: 

ansatz density profile n(r) = no(l — e r 1 ) where A is picked 
so that the two densities agree at one coherence length away 
from the vortex center. 

PHY-0456720 and PHY-0803371, and The Research Cor- 
poration Cottrell Scholars program (GR). 

APPENDIX: DENSITY PROFILE FOR A SINGLE 
VORTEX 

In order to find the short-ranged interspecies vortex 
interaction, we need to know the behavior of the den- 
sity of the condensate about a vortex. To this end we 
consider the Gross-Pitaevskii equation for a single com- 
ponent BEC having a vortex at the origin. That is, we 
write ip = fe %0 and take 6 = cp where ip is the azimuthal 
angle from polar coordinates. Substituting this into the 
GPE leads to the following equation dictating the density 
profile 

-£>W) + |^ + 5 /3 = ,/. (A.1) 



where A is a parameter on the order of the coherence 
length. Note that while this ansatz has the correct form 
close to the vortex core, the long distance behavior differs 
considerably. Fortunately our problem of vortex lock- 
ing is dominated by the short-distance behavior, and we 
choose A so that the two densities (numerical and ansatz) 
agree at r = £ which requires A = 1.781£, as shown in 
Fig. El 

Our positive results, confirming the vortex-locked 
state, also confirm our intuition that the origin of the 
phenomena is in the short-ranged attraction between vor- 
tices. As explained in Ref. Q] (but without proof), the al- 
gebraic decay of the superfluid order parameter of a single 
votex does not imply that vortices of one species, when 
in a lattice, exhibit a power-law decaying force on the 
vortices on the other species. Unlike the single-species 
vortex- vortex force, which is the result of the inductive 
(kinetic) energy term in the Gross-Pitaevskii equations, 
the interspecies force is a result of a density-density in- 
teraction. The density suppression due to a single vortex 
occurs since the superflow of the vortex effectively in- 
creases the mass terms Vi, V2 in Eq. (jlj). But in a lattice 
of vortices, the combined superflow vector is nearly zero 
(i.e., negligible compared to 7i/ra a £ a ), and, therefore, so 
is the respective density suppression. 
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